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q>^ Abstract 

£S) We measure the scaling properties of the probability distribution of the smoothed den- 

sity held in .AT-body simulations of expanding universes with scale-free initial power-spectra, 
(l^fc| 2 ) x k n , with particular attention to the predictions of the stable clustering hypothesis. 
nI We concentrate our analysis on the ratios Sq(1) = £q/£,2 where £q is the averaged Q-body 

correlation function over a cell of radius I. According to the stable clustering hypothesis, Sq 
should not depend on scale. We measure directly the functions Sq(£) for Q < 5. The behavior 
of the higher order correlations is studied through that of the void probability distribution, Po, 
ON which is the probability of finding an empty cell of radius I. If the stable clustering hypothesis 

applies, the function Po should also exhibit remarkable scaling properties. 
O | In our analysis, we carefully account for various spurious effects, such as initial grid con- 

' tamination, loss of dynamics due to the short range softening of the forces, and hnite volume 

; I size of our simulations. Only after correcting for the latter do we hnd agreement of the mea- 

sured Sq, 3 < Q < 5 with the expected self-similar solution Sq(i,t) = Sq(£ 2 ) = Sq{1/ lo{t)), 
C$ io(t) <x i 4 /( 9 + 3n ). The void probability is only weakly sensitive to such defects and closely 

follows the expected self-similar behavior. 

As functions of £ 2 , the quantities Sq, 3 < Q < 5, exhibit two plateaus separated by a smooth 
transition around £ 2 ~ 1- I n the weakly nonlinear regime, £ 2 1) the results are in reasonable 
agreement with the predictions of perturbation theory. In the nonlinear regime, £ 2 > 1) the 
function Sq(£ 2 ) is larger than in the weakly nonlinear regime, and increasingly so with —n. It 
is well-htted by the expression Sq = (^ 2 /100) ' 045 ^ -2 ' Sq for all n. This weak dependence on 
scale proves a small, but significant departure from the stable clustering predictions at least for 
n = and n = +1. It is thus also at variance with the predictions of the hierarchical model. 

The analysis of Po confirms that the expected scale-invariance of the functions Sq is not 
exactly attained in the part of the nonlinear regime we probe, except possibly for n = —2 and 
marginally for n = — 1. In these two cases, our measurements are not accurate enough to be 
discriminant. On the other hand, we could demonstrate that the observed power-law behavior 
of Sq cannot be generalized as such to arbitrary order in Q. Indeed this would induce scaling 
properties of Po incompatible with those measured. 

subject headings: galaxies: clustering - methods: numerical - methods: statistical 
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1 Introduction 



Large-scale structures in the observed galaxy distribution are thought to have arisen from small 
initial fluctuations through gravitational instability. Generally, it is assumed that the large-scale 
dynamics is dominated by collisionless dark matter. Of course, galaxy formation is not a collisionless 
process, and the extent to which one can infer some properties of the distribution of luminous matter 
from that of dark matter is a matter of debate. Moreover, statistical indicators in observational 
catalogs are subject to many contamination effects, such as the finiteness of the sampled volume, 
selection effects, redshift distortion in three-dimensional samples, and projection effects along the 
line of sight in two-dimensional sky surveys. 

Albeit simpler than galaxy formation, the dynamics of collisionless matter on large scales in 
the Universe is still not fully-understood. The statistical evolution of self-gravitating collisionless 
particles may be described by the BBGKY hierarchy (e.g., Peebles 1980, hereafter LSS), which is 
an infinite system of coupled differential equations for the (J-body correlations functions in phase- 
space. In the weakly nonlinear regime, where the fluctuations of the density field are small, the 
hierarchy can be closed through perturbative approaches. However, it does not appear possible 
to solve the BBGKY hierarchy analytically when fluctuations of the density field become large, 
owing to the long-range nature of gravitational forces. In the highly nonlinear regime, therefore, 
JV-body simulations are required to study the evolution of the density distribution and to measure 
the statistical properties of the system under consideration. 

In this paper we analyze flat universes consisting entirely of collisionless matter and having 
scale-free initial Gaussian fluctuations 

P(k) = (\6 k \ 2 ) = Ak n . (1) 

To analyze the statistics of gravitational clustering, we employ the count probability distribution 
function P^{m,l) (CPDF), which represents the probability of finding N objects in a cell of radius 
I located randomly in a sample of points with average number density m. The generating function 
V{\) = J2^ N Pn of the CPDF can be written (White 1979; Schaeffer 1984; Balian & Schaeffer 
1989a, hereafter BS; Szapudi & Szalay 1993) 



V(\) = exp < 
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where £q is the average of the (J-body correlation function £q (e.g., LSS) over a cell: 
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The CPDF and its moments are thus well-suited to describing the scaling properties of a system 
of particles. From the dynamical point of view, an especially interesting quantity is the ratio 

SQW = (4) 
£2 

Indeed, as discussed below, Sq is expected to vary only weakly with scale or even to obey the 
following scaling relation over some range 

Sq = constant of scale, (5) 
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as suggested by measurements of the low-order correlation functions £q in the observed galaxy dis- 
tribution (e.g., Groth & Peebles 1977; Fry & Peebles 1978; Davis & Peebles 1983; Sharp, Bonometto 
& Lucchin 1984; Szapudi, Szalay & Boschan 1992, Bouchet et al. 1993). Note that S3 and S4 are 
merely the (renormalized) skewness and kurtosis of the smoothed density distribution function. 

In principle, it is straightforward to measure Sq because this quantity is simply related to 
the moments of the CPDF (e.g., Szapudi & Szalay 1993). However, for larger and larger Q, Sq 
is increasingly dominated by the high N tail of the CPDF which is subject to spurious sampling 
effects related to the finite size of the volume (Colombi, Bouchet & Schaeffer 1994, hereafter CBSI). 
It is possible to account for these difficulties and accurately determine Sq, but only for moderate 
orders Q & 10 (CBSI). In what follows, we limit our analysis to the fifth order, Q = 5. Higher- 
order statistics are considered indirectly through the void probability distribution function Po(m,l) 
(VPDF). From equation (2), it is indeed easy to write Po(m,l) as 

P {m,l) = exp [-JV<T(m,£)] , (6) 

where 

aK£)=E(-l) Ar - 1 ^ 1 ^- 1 - (7) 
N>1 

The quantity W = (4/3)7rm£ 3 is the average number of objects per cell, and JV C is the typical 
number of objects in a cell located in an overdense region 

N C = N^ 2 . (8) 

It is natural to study the scaling behavior of the VPDF by determining a as a function of JV C . 
Indeed, if the scaling relation (5) applies, then (White 1979, BS) 

a(n,l) = a(N c ). (9) 

Measurements of the VPDF in the observed galaxy distribution are in good agreement with equa- 
tion (9) (e.g., Sharp 1981; Bouchet & Lachieze-Rey 1986; Maurogordato & Lachieze-Rey 1986; Fry 
et al. 1989; Maurogordato, Schaeffer & da Costa 1992). 

In the weakly nonlinear regime (WR), £ 2 "C 1, the equations of motion can be solved perturba- 
tively. For example, second-order perturbation theory (LSS) yields £3 (Fry 1984a) and hence the 
skewness 5s(£) of the smoothed density field (Juszkiewicz & Bouchet 1992; Juszkiewicz, Bouchet & 
Colombi 1993, hereafter JBC). Higher-order perturbation theory can be used to compute additional 
correlations (Goroff et al. 1986; Bernardeau 1992), such as the kurtosis S^l) (Bernardeau 1994a; 
Lokas et al. 1994). In fact, the full hierarchy of averaged correlations £q(1) can be determined 
in the limit £ 2 "C 1 (Bernardeau 1994b, hereafter B94). Note that the hierarchy of the ratios Sq 
depends only weakly on the cosmological parameters and A, at least for Q = 3 (Bouchet et al. 
1992, 1995; Hivon et al. 1995). As shown below, for scale-free initial conditions and = 1, the 
function £ 2 1S a power-law of scale and Sq obeys the scaling relation (e.g., JBC, B94). 

In the strongly nonlinear regime (SR), £ 2 ^ 1) no general analytical solutions to the BBGKY 
hierarchy have been found. Indeed, it is hard to deal with an a priori infinite hierarchy of correla- 
tions, since in this regime £ 2 <C £3 <C ••• <C £q, so standard truncation methods cannot be applied. 
However, there is no preferred scale in gravitational dynamics and the BBGKY hierarchy admits 
self-similar solutions (e.g., Davis & Peebles 1977). Moreover, local statistical equilibrium should 
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obtain on small enough scales, corresponding to virialized objects. In that case, the BBGKY hier- 
archy simplifies and the ratios Sq obey the scaling relation (e.g., Davis & Peebles 1977; Balian & 
Schaeffer 1989b). However, except in some particular frameworks (Fry 1984b, Hamilton 1988), the 
values of Sq in this regime are currently unknown and only measurements in JV-body simulations 
can help to determine them. In particular, there is no reason why the values of Sq should be the 
same in the SR and in the WR. 

In existing three-dimensional galaxy catalogs, the difference between these two regimes is weak, 
i.e. the property (5) seems to apply over all the available dynamic range (e.g., Bouchet et al. 1993), 
but this is likely to be due, at least partly, to projection in redshift space which tends to flatten the 
functions Sq(1) (Lahav et al. 1993; Matsubara & Suto 1994; Hivon et al. 1995). A bias between 
the distribution of galaxies and the mass could also flatten the Sq. In two-dimensional galaxy 
catalogs, the Sq show a significant scale dependence (Gaztahaga 1994), but the interpretation of 
this result is complicated by projections effects. 

Measurements of low-order correlations in JV-body simulations indicate that Sq changes from 
the WR to the SR (e.g., Bouchet & Hernquist 1992, hereafter BH; CBSI; Lucchin et al. 1994). 
However, while the WR agreement between JV-body simulations and perturbation theory is well- 
established (e.g., JBC; B94; Gaztahaga & Baugh 1995), the nonlinear regime behavior is still 
quite uncertain. The early JV-body experiments of Efstathiou et al. (1988, hereafter EFWD) 
with scale-free initial conditions indicate that the (J-body correlations (Q < 3) exhibit self-similar 
and scale-invariant behavior in the SR, in agreement with Fry, Melott & Shandarin (1993). In 
a similar spirit, if the scaling relation applies, the CPDF itself should display remarkable scaling 
properties, as predicted by BS. Although testing these predictions is quite a delicate matter, the 
scaling relation can be strongly discriminated against in this way (CBSII). Bouchet, Schaeffer & 
Davis (1991, BSD) and BH measured the CPDF in JV-body simulations with Cold Dark Matter 
(CDM), Hot Dark Matter (HDM) and white noise initial conditions and found striking agreement 
with the predictions of BS. CBSI measured the functions Sq(1), Q < 5 in the same simulations and 
found good agreement with scale-invariance in the highly nonlinear regime, but only after correcting 
for finite volume effects. However, there is not yet a consensus in the scientific community as to 
whether or not the scaling relation is reached in the highly nonlinear regime. For example, Lahav 
et al. (1993) and Suto & Matsubara (1994) who measured respectively the low-order moments of 
the CPDF and the functions £q, Q < 4 in low density CDM and white noise simulations, find a 
significant disagreement with equation (5) in the SR. 

The measurement of correlations is difficult because the samples we have access to are not 
perfect. First, the finite volumes of the samples can induce unphysical distortions of the Sq 
(CBSI). Other spurious effects, such as discreteness, or artificial correlations related to a particular 
way of setting up initial conditions (CBSII), as well as other numerical effects due to some intrinsic 
properties of the JV-body codes used to make the simulations can also introduce systematic errors. 
One has to carefully control all these defects before accepting any far-reaching conclusions. 

In this paper, we carefully measure the quantities Sq(1) for Q < 5 and the VPDF in JV-body 
simulations with scale-free initial power-spectra given by equation (1) with n = —2, —1, 0, +1. We 
explicitly account for the possible contamination effects noted above, particularly those related to 
the finite volume of the simulations. In addition to exploring the questions of whether or not the 
scaling-relation is obeyed in the highly nonlinear regime, our simulations enable us to study how 
the scaling behavior of the system changes with initial conditions. 

The paper is organized as follows. In § 2, we describe the JV-body simulations and how we 
measured the CPDF in our samples. The self-similar solution expected in our scale-free systems is 
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recalled in § 3, along with possible spurious effects which can induce artificial deviations from such a 
solution. In § 4 we measure the two-body correlation function, compare the results with theoretical 
predictions, and determine a scale range where the effects noted in § 3 should be negligible or can 
be accounted for. We then measure the ratios Sq, Q < 5 and show the necessity of correcting for 
finite volume effects. We show that agreement with perturbation theory is good in the WR. In the 
SR, the expected scale invariance does not seem to be exactly achieved, except perhaps for n = — 1 
and n = —2. The deviations from equation (5) are however very weak. Section § 5 deals with the 
void probability distribution. We test the agreement of the measured VPDF with the self-similar 
behavior and discuss all the possible contamination effects introduced in § 3. We can then carefully 
measure a as a function of JV C in the nonlinear regime and see if our conclusions for the low-order 
correlations can be generalized to higher order. In § 6, we summarize the results and conclude the 
paper with a short discussion. 

2 The iV-body samples 

We performed a set of five scale-invariant simulations of flat universes involving N paT = 64 3 dark 
matter particles in a cubic box of width £t>ox = 1 with periodic boundaries using the cosmological 
treecode of Hernquist, Bouchet & Suto (1991). The tolerance angle was 9 = 0.75 and the softening 
parameter e of the short range component of the force was given by e = A par /20, where 

Apar = L hoK N-iJ 3 ~ 0.016 (10) 

is the mean interparticle distance. The initial power-spectra are given by equation (1) with n = —2, 
— 1, and +1. Two simulations with different initial random phases were performed for the case 
n = —1. Following EFWD, we chose the normalization factor A so that the amplitude of the 
power-spectrum matched the white-noise level at the Nyquist frequency of the particles, except for 
n = —2. In that case, A was taken four times smaller. The timestep choice is discussed in Appendix 
A. Table 1 lists the values of the expansion factor for which we have analyzed the JV-body data. 

For each simulation and each expansion factor analyzed, we measured the CPDF for spherical 
cells of radius I in the following scale range (in units of £t>ox) 

- 2.8 < \og w l < -1.0, (11) 

with a logarithmic step, Alog 10 £ = 0.2. The lower limit l\ OVf is imposed because we want to sample 
scales larger than the spatial resolution of the simulations, which should be a few e. We took 
l\ OVf ~ 2e, but we shall see later that this constraint is not sufficiently strong. The upper scale 
4nax was chosen so that the number of independent modes in Fourier space corresponding to l max 
is large. This does not guarantee, however, that finite volume effects are negligible, as we shall see 
later. 

To compute the CPDF, we sample our simulations with a regular pattern of cells, involving 
JV no des(-0 points. Table 2 gives the chosen iV no d es f° r all the snapshots except for (n, a) = ( — 2, 2), as 
a function of scale and the corresponding minimal value of the CPDF P™ n that can be measured. 
It is important to realize that, to have an accurate sampling, the mean inter-cell separation must 
be much smaller than the scale of interest, which is the case with the choice of iV no d es given in Table 
2. The statistical errors arising from the finite number of cells used to sample the JV-body catalogs 
are quite small, in fact negligible compared with the other errors, such as those related to finite 
volume effects. For the case (n, a) = ( — 2,2), we used a smaller number of cells than for the other 
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Table 1: Important scales in the simulations 



n a 
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s d 




logio4 


lo SlO^M 


FV h 


symbol 1 




2.0 


0.021 


0.021 


0.0048 


TT 

-2.0 # 


-1.0 


N 


square 


-2 


3.2 


0.051 


0.052 


0.0021 


_2.4#> @ 


-1.4* 


Y 


triangle 




5.2 


0.103 


0.138 


n.a. 


-2.4® 


-1.0 


Y 


star 




8.0 


0.173 


0.327 


n.a. 


-2.4® 


-1.2' 


Y 


diamond 




2.5 


0.022 


0.022 


0.0029 


n 

-2.2* 


-1.0 


N 


square 


-1 


6.4 


0.057 


0.054 


n.a. 


-2.4® 


-1.4* 


Y 


triangle 




16 


0.131 


0.136 


n.a. 


-2.4® 


-1.0 


Y 


star 




4.0 


0.019 


0.019 


0.0022 


-2.4#>® 


-1.0 


N 


square 





16 


0.048 


0.047 


n.a. 


-2.4® 


-1.4* 


Y 


triangle 




64 


0.118 


0.120 


n.a. 


-2.4® 


-1.0 


Y 


star 




6.4 


0.020 


0.020 


0.0016 


-2A* 


-1.0 


N 


square 


+ 1 


41 


0.048 


0.049 


n.a. 


-2.6® 


-1.4 


N 


triangle 




256 


0.119 


0.123 


n.a. 


-2.6® 


-1.0 


Y 


star 



a spectral index. 

b expansion factor (assuming that initially a = 1). 
c the measured correlation length. 

d the correlation length according to the self-similar solution (scaled to Iq at a = 2.0, 2.5, 4.0 and 
6.4 respectively for n = —2, —1, and +1). 

e the typical distance between two matter particles in an overdense region, when measurable. 

f minimum reliable scale. The indices indicate the effects imposing I > l m , * for grid effects (we 
impose I > 1.5£ c ) and ® for finite force resolution effects. 

g maximum available reliable scale, when measuring and correcting for finite volume effects for the 
quantities Sq, Q = 3, 4, 5. The indices indicate when some constraint has to be obeyed, * 
means that I > Iq was imposed for the correction to be valid, 7r means that larger scales were 
two noisy to be able to meaningfully correct the quantities Sq for finite volume effects. 

11 indicates if finite volume effects (FV) have been corrected for (Y) or not (N) in figure 4. 

1 symbols used in the figures. 

cases, iV no des = 512 3 for log w i < —1.4 (so we probably slightly undersample for log w i = —2.8), 
iV nod e S = 384 3 for log w l = -1.2 and JV nodes = 128 3 for log w l = -1.0. 

3 Self-similarity and numerical limitations 

Since there are no preferred scales in the initial power-spectra, there should be only one relevant 
physical scale in the simulated system: the correlation length Iq defined by 

UIq) = 1. (12) 

This scale separates the weakly nonlinear regime (WR) where £2^1 fr° m the strongly nonlinear 
regime (SR) where £ 2 ^ 1- -A- n Y statistical quantity f(i,t) can thus be written as 

f(l,t) = g(l/l ). (13) 
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Table 2: Resolution of the sampling cells pattern as a function of logarithm of scale and the 
corresponding minimal CPDF reached 



logio^ 


-2.8 


-2.6 


-2.4 


-2.2 


-2.0 


-1.8 


-1.6 


-1.4 


-1.2 


-1.0 


nodes 


2176 


1792 


1536 


1280 


1024 


896 


768 


640 


512 


384 


lncr pmin 


-10.0 


-9.75 


-9.55 


-9.3 


-9.0 


-8.85 


-8.65 


-8.4 


-8.1 


-7.75 



In the WE,, the two body correlation function is £ 2 = C(n)a?l n 3 , where C is a constant depending 
only on the initial power spectrum index (see § 4.1.1 for the exact expression of £ 2 )j which implies 
that (Inagaki 1976, Davis & Peebles 1977, LSS) 

to = s(t), s(t) oc a a , (14) 

with 



The third column of Table 1 gives the measured correlation length in units of £boxj to be compared 
to the values of s(t). The latter are given in the fourth column of Table 1, assuming that the initial 
value of s(t) is exactly equal to the measured Iq at the smallest sampled expansion factor. The 
agreement of the measurements with equation (14) is very good, except for the case n = —2, where, 
at late stages, the deviation from self-similar behavior appears at first sight quite significant (for 
a <; 5.2). We shall see in § 4.1.3 that this lack of agreement with self- similarity is spurious and is 
due to finite volume effects. 

Deviations from the expected behavior given by equations (13) and (14) should only be due to 
numerical limitations of JV-body simulations, unless linear theory is not valid in the WR. Here, as 
in Hivon et al. (1994), we list numerical limitations of JV-body simulations that can contaminate 
the measurements. In § 3.1, we discuss possible dynamical effects due to the discreteness of our 
representation of an underlying continuous density field. § 3.2 deals with grid effects and transients 
resulting from the use of the Zel'dovich approximation to slightly perturb a regular pattern of 
particles to set up initial conditions. In § 3.3 we mention a possible lack of accuracy at small scales 
due to the short range softening of the forces in the simulations. In § 3.4, we discuss finite volume 
effects. 



3.1 Collisionless fluid limit 

The JV-body code should describe the evolution of a "self-gravitating" system of "micro" -particles. 
On the scales of interest, the collisionless (mean-field limit) Boltzman equation, i.e. the Vlasov 
equation, is valid (see, e.g., LSS). The "macro" -particles we use in JV-body codes have masses 
much larger than those of the micro-particles. Even if the real underlying distribution is well 
approximated by a mean-field description, it may not be the case for our system of macro-particles 
in the scaling regime we consider. To reduce dynamical effects arising from the discreteness of 
our representation, such as two-body relaxation, a short range softening parameter e is introduced 
(see § 3.3, Appendix A). In itself, however, this is not sufficient. Before making measurements 
on the JV-body simulation, one must wait long enough for the numerical system to indeed achieve 
the mean-field limit of a "collisionless fluid". In other words, the typical number of particles per 
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collapsed object (of size similar to the Jeans length) must be large. The typical size of a clump is 
the correlation length. Thus, we should have 

to > A par ~ 0.016. (16) 

This condition is fulfilled for the expansion factors we chose to analyze, although barely for the 
smallest ones (corresponding to a = 2.0, 2.5, 4.0 and 6.4 for n = —2, —1, and 1 respectively, see 
Table 1). 

3.2 Initial conditions 

To reduce the small-scale shot-noise arising from the discreteness of our "macro" -particles, one 
traditionally starts the simulations from a regular distribution of particles, slightly perturbed by 
using the Zel'dovich approximation (Zel'dovich 1970) to displace the particles. Some imprint of this 
initial pattern is conserved during the simulation, particularly in underdense regions. It can thus 
contaminate the measurement of the CPDF, especially at small N (e.g., BSD, BH, CBSII). From 
the statistical point of view, the discrete realization of a continuous density field should indeed be 
locally Poissonian. We shall discuss in more detail such "grid" effects on the VPDF in § 5. As far as 
the low-order correlation functions are concerned, these effects should be negligible, except in the 
early stages of the simulation, at scales smaller than, or of the same order as the mean interparticle 
distance l c in overdense regions. 

The low-order correlations are dominated by regions with high density contrasts, and increas- 
ingly so with Q. Because nonlinear gravitational dynamics is intrinsically chaotic, it is likely that 
such regions become locally Poissonian after shell-crossing. To have such phase-space mixing, con- 
dition (16) must be fulfilled, but it is not obvious that this is a sufficient condition: it is possible to 
have high density contrasts without shell-crossing, particularly if there is a cut-off at small scales 
in the spectrum of initial fluctuations. Note that since the ratio between small-scale power and 
large-scale power decreases with —n, we expect grid effects to be less important for large values of 
n. They should also become less and less significant as the system evolves. 

When the typical number of particles per cell located in an overdense region is large compared 
to unity, grid effects should be negligible. This condition yields a constraint on the scales which 
are fully reliable of the form 

l>l c . (17) 
The scale l c is the typical distance between two particles in a cluster, defined by (BS) 

N c {l c ) = 1. (18) 

Some measured values of l c , when available, are given in Table 1. We transformed the asymptotic 
constraint (17) in the more practical condition 

I > 1.5l c . (19) 

We see for example that, for the first snapshot, we should have log 10 l > —2.6, —2.50, —2.35 and 
— 2.15 respectively for n = 1, 0, —1, —2 for grid effects to have a negligible influence on the low-order 
correlations. This partly explains our choice of the "minimum reliable scale" l m given in Table 1. 
For larger expansion factors, the measured clustering number JV C is always sufficiently large in the 
available dynamic range (11) for grid effects to be insignificant (except for (n, a) = ( — 2,3.2) where 
we should also have log 10 l > —2.5). 
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Figure 1: The CPDF measured in the most evolved snapshot of our simulations as a function 
of N, for various scales I. The system of coordinates has been chosen to emphasize the large- N 
tail of the CPDF, which is clearly exponential, although rather noisy. As N increases, the CPDF 
indeed presents larger and larger irregularities followed by a sharp cutoff. The smooth lines are 
the analytical fits we derive in § 4.2.2. In each panel, the top curve corresponds to the scale 
\og w l = —1.0. Going downwards, the scale decreases with each curve, with a logarithmic step 
Alog 10 £ = 0.2. 

The use of the Zel'dovich approximation to set up initial conditions can also induce some spu- 
rious effects (e.g., JBC). Indeed, this approximation describes a first (linear) order approximation 
to the trajectories of the particles. It is not valid if one aims at measuring higher order quantities, 
such as the skewness S3 of the density distribution. It is thus necessary to wait sufficiently long that 
higher-order coupling terms had time to fully develop. The corresponding requirement is a> 1, 
which is hardly fulfilled by our simulations. Practically, however, a <; 3 should be sufficient (e.g., 
Baugh, Gaztahaga, & Efstathiou, 1994). We thus expect the first snapshot of our simulations to be 
contaminated by transients, particularly on large, weakly nonlinear scales. Note that the function 
£ 2) which is of first order in the WR, should be only weakly influenced by these effects. 

3.3 Short range softening of the forces 

The dynamics on the smallest scales in our simulations is not accurate because of the softening 
parameter e used to bound forces and reduce two-body relaxation (§ 3.1). Typically, however, if the 
timestep is carefully chosen (see, e.g., Appendix A), the softening parameter should contaminate 
measurements up to at most a few e. We shall see, however, that softening can affect the count-in- 
cells statistics on scales as large as \og w l ~ —2.5 (§ 4.1.2), which partly explains our choice of the 
"minimal reliable scale" l m in Table 1. 

3.4 Finite volume effects 

Since we use periodic boundaries in our simulations, the fluctuations from scales larger than the 
box size are missing. Also, on scales smaller than but comparable to the box size, only a few 
independent modes of the power spectrum are sampled, thereby increasing the uncertainties in 
measurements as one approaches the box size. To minimize these effects, one usually requires 
log 10 l Sa —1, as we did. But this constraint is not necessarily sufficient. Indeed, the large N tail 
of the CPDF is determined by just a few large clusters and is thus subject to fluctuations due to 
small number statistics, until it reaches an artificial cutoff at iV max (CBSI, CBSII), as illustrated 
by Figure 1. If the sample size £t>ox is large enough compared to the correlation length (typically 
-^box *S 20^o)) which is the typical size of a cluster, these modifications of the true tail occur for 
such small values of the CPDF that they are of no consequence. When the correlation length 
becomes significant compared to £t>ox (typically Iq <; Zbox/20), this scale-dependent effect is likely 
to influence measurements of the low-order correlation functions £q, increasingly with Q. Also, 
this finite volume effect should be more significant when the amplitude of the fluctuations on large 
scales is big, so it should increase with —n. Fortunately, it can be corrected for, or at least an 
error-bar associated with it can be evaluated (CBSI). 
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4 Scaling behavior of the low-order correlations 



4.1 The averaged two-body correlation function 
4.1.1 Theoretical predictions 

In the WE,, the function £ 2 can be written (e.g., Hivon et al. 1994) 



q on-i r(l -n)T P±s, 

Ui) = — \ 2 ' Aa 2 r n - 3 , I > to (20) 

521 1 4tt T 2 ^ r ^ 1 1 

where A is the initial amplitude of the power-spectrum and a is the expansion factor, assumed to 
be unity at the beginning of the simulation. 

In the SR, with the additional hypothesis of local statistical equilibrium, the two-body correla- 
tion function should be a power-law of scale (Davis & Peebles 1977): 

£ 2 (f)ocr 7 , Kl , (21) 

with 



7 = (22) 

' 5+n K 1 

4.1.2 Measurement 

Figure 2 displays the quantity £ 2 as a function oil/ s{t) (where s(t) is given by eq. [14]) in logarithmic 
coordinates and for the various expansion factors listed in Table 1. To compute £ 2 we use the 
following formula (Fry & Peebles 1978, LSS), that corrects for discreteness 

N%{1) = H2-N, (23) 

where /x 2 is the centered, second moment of the CPDF. More generally, the centered moment of 
order Q is defined by 

oo 

^(1) = {(N -Nf)=J2(N- N) Q PnW- (24) 

N=0 

We show results for only one of the two n = — 1 simulations. Indeed, the measured £ 2 are nearly 
identical in the two simulations, except perhaps for the largest expansion factor a = 16, for which 
the differences are anyway smaller than A£ 2 /£ 2 ^ 0.3. This number corresponds to a vertical 
length in figure 2 approximately 1.5 times smaller than the size of the symbols used to make the 
plots. 

As displayed in Figure 2, if the evolution of the system is self-similar, all the curves should 
superimpose. Possible deviations from such behavior must be induced by numerical effects. The 
short dashed lines indicate the logarithmic slope expected from equations (21) and (22) in the SR. 
The normalization has been chosen to match the measurements. The long dashes give the linear 
theory prediction (20), expected to apply in the WR. 

Except for the case n = —2, which we discuss further, the measured £ 2 (^,i) is in very good 
agreement with theoretical predictions. Deviations from self- similarity are only significant on the 
smallest scales, typically for log 10 l ^ l m , with 

l m ~ -2.5, (25) 
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Figure 2: The averaged two-body correlation function £ 2 as a function of l/s (where s is given 
by eq. [14]) in logarithmic coordinates, measured in each simulation for various expansion factors 
(Table 1; see also this table for the significance of the symbols). The long dashes give the linear 
prediction, valid in the limit £2 ^ 1) an d the short dashes indicate the logarithmic slope of £ 2 
expected in the limit £ 2 ^ 1 (eqs. [21] and [22]). In the last case, the normalization has been 
chosen to fit the data. The dotted curve in the bottom right panel is a polynomial fit in logarithmic 
coordinates (see eqs. [27] and [28]). 
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which corresponds to the two left-most points of each curve in each panel of Figure 2, for which 
the measured £ 2 underestimates the true value. This effect seems to increase with — n and tends to 
decrease at late stages in the simulations. For the first expansion factor analyzed in each simulation 
(or the first two for n = —2), the deviation is mainly due to the numerical limitations discussed 
in § 3.1, 3.2 and 3.3, particularly grid effects. For larger expansion factors, grid effects become 
negligible and only the short range softening of the forces introduced in § 3.3 should matter. 

A careful measurement of the two-body correlation function ^( r ) shows that at late stages 
in the simulations the system exhibits self-similar behavior on scales as small as l e ~ 2e. This 
excellent agreement with theoretical predictions does not apply for velocity correlations, which are 
much more sensitive to the force resolution. Note that the function £ 2 (-0 1S related to ^( r ) through 
an integral over a cell (eq. [3]). In this integral, small separations |ri — r 2 | 5a l e can contribute 
significantly. Only when I is large enough compared to l E , will the effect of softening become 
negligible. Experimentally, we find in equation (25) that £ 2 (-0 can be contaminated up to the scale 
l m ~ 2l s ~ 4e. Table 1 gives more accurate estimates of the "minimum reliable scale" l m we chose. 

4.1.3 Finite volume effects 

The very good agreement with a self-similar behavior on scales larger than l m , and the convergence 
to the predictions of linear theory on large scales 1 indirectly suggest that finite volume effects on £ 2 
are negligible in the scaling range we consider, for n = 1, and —1. This can also be checked in a 
different manner with the two-body correlation: as discussed in CBSI, because we miss some power 
on scales larger than the simulation box size, the measurements of £ 2 (-£) are likely to underestimate 
its true value, increasingly with I. In a box of finite volume, the measured function ^(-O can be 
written 

— v, 7"^ f°° 

^2° X W~^f J 1/3 ^t° x \ 2 )[We(k)] 2 k 2 dk, (26) 

where (|#J? ox | 2 ) is the power spectrum measured in the box and Wi = 3(kl)~ 3 [s'm(kl) — klcos(kl)] 
is the Fourier transform of the top hat filter. One can correct such measurements for finite volume 
effects by extrapolating the measured power-spectrum (|#J? ox | 2 ) on scales larger than £t>ox using 
linear perturbation theory and by computing the integral (26) in the limit £t>ox —> 00 • Of course, 
such a correction is valid as long as the large-scale dynamics does not dramatically influence the 
power-spectrum on scales smaller than the box, so the correlation length should be a reasonably 
small fraction of the box size. We applied the above procedure to our data and indeed noticed that 
it yields only small changes for n > —2. The correction is significant for n = — 2 for all a when 
log 10 l > —1.4, but does not increase £ 2 by more than a factor of order of 10 ' 1 at log 10 l = —1.0. 
This result is in agreement with the bottom right panel of figure 2 for a = 2 and a = 3.2. But for 
a = 5.2 and a = 8.0, it appears that finite volume effects are much stronger than our correction 
seems to indicate. However, in that case, our prescription is unreliable, because the correlation 
length becomes a significant fraction of the box size, larger than Zbox/10, so nonlinear coupling can 
affect Fourier modes with wavelengths as large as the size of the box or even larger, particularly at 
a = 8. 

The strong deviation from self- similarity on scales larger than l m in the n = —2 case is thus due 
to finite volume effects. The true function £ 2 is larger than the measured one, increasingly with 
scale. The correction invoked above shows however that one can reasonably infer from Figure 2 

1 except for n = +1: in this case, nonlinear effects might never be negligible, whatever scales are considered. 
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the global shape of £ 2 , assuming that at each scale the true value of £ 2 1S approximated (within 
a factor of 10 ' 1 ) by the maximum of the values obtained in each of the snapshots analyzed. The 
corresponding dotted curve in Figure 2 is given by the following function F(x), in logarithmic 
coordinates 

F(x) = ai(cc - x{) + 2/1, x < x 1: 

F{x) = P 3 (x), Xl <x< x 2 , (27) 

F{x) = a 2 (x - x 2 ) + 2/2, x > x 2 , 

where Pz(x) is the polynomial of degree four satisfying 

P3(x 1 ) = y 1 , dP 3 /dx 1 = ai, P 3 (x 2 ) = y 2 , dP 3 /dx 2 = a 2 , P 3 (x 3 ) = y 3 , (28) 
with xi = -1.7, j/i = 2.42, a x = -1, x 2 = 0.6, y 2 = -7.4, a 2 = -1, x 3 = -0.7, y 3 = 1.11. 



4.2 Higher order correlation functions 
4.2.1 Theoretical predictions 



The measurement of £ 2 allows us to define a range of scales for which the JV-body simulations for 
n > —2 exhibit the expected self-similar behavior given by equations (13) and (14). Furthermore, 
in the SR, £ 2 is compatible with the predictions obtained when local statistical equilibrium is 
assumed. If this latter hypothesis is valid, the quantities Sq should obey the scaling property (5) 
in the regime £ 2 3> 1. In the WR regime, the quantities Sq are expected to match the predictions 
of perturbation theory, i.e., 

34 

S 3 = — -(n + 3), (29) 
(Juszkiewicz & Bouchet 1992, JBC), and (B94) 

60712 62, oN 1, oN2 
Si = 7323" " Y {n + )+ 3 {n + 3) ' (30) 

200575880 1847200, . 6940, x2 235, x , , . 

S 5 = (n + 3) + (n + 3) 2 (n + 3) 3 . (31) 

305613 3969 K J 63 K J 27 K ' K ' 

4.2.2 Measurements: correcting for finite volume effects 

Figure 3 displays the ratios Sq = £q/£ 2 1 as functions of £ 2 as measured in our simulations, in 
the available scale range given by equation (11). To compute the averaged correlation functions 
£q, we use the standard formulae, that correct for discreteness (Fry & Peebles 1978, LSS, CBSI) 
and generalize to higher orders from equation (23): 

= Ms - 3/x 2 + 2N, (32) 

~N%(1) = m- 6/x 3 - 3fi 2 2 + Un2 ~ 6N, (33) 
= Ms - 10/x 4 - 10/X2M3 + 35/x 3 + 30fi 2 2 - 50/x 2 + 24JV. (34) 

The quantity [Lq is the centered moment of order Q of the CPDF (eq. [24]). 

Figure 3 shows that the Sq do not obey self- similarity as well as the averaged two-body corre- 
lation function, and it becomes worse with decreasing n. All the curves corresponding to various 
expansion factors a but the same value of Q (at fixed n) should superimpose, which is far from 
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Figure 3: The quantities Sq = £q/£ 2 as functions of £ 2) i n logarithmic coordinates, measured 
in each simulation for various expansion factors (Table 1; see also this table for the meaning of the 
symbols) in the full available scaling range (eq. [11]). In each panel, the values of Sq increase with 
the order Q. In this coordinate system, since there is no preferred scale, all the curves (at fixed Q 
and n) should superimpose, which is far from being the case, even more so as — n increases. 
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being the case, particularly for n = —2. At a fixed value of £ 2 , the measured value of Sq decreases 
with time. These effects are spurious, due to the finiteness of the sampled volume. As discussed 
by CBSI and in § 3.4, finite volume effects are more important when the correlation length is a 
large fraction of the sample size and when there is more power at large scales. They induce in- 
creasing irregularities with N on the CPDF until a sharp cutoff at N = N maK (l) is encountered, 
as illustrated by Figure 1. These irregularities, and particularly the cutoff at iV max contaminate 
the measurement of the low-order averaged correlation functions, which are increasingly sensitive 
to the large- N part of the CPDF as Q increases. 

Fortunately, it is possible to use the procedure proposed by CBSI to correct for finite volume 
effects, by smoothing and extending to infinity the large- N exponential tail exhibited by the CPDF. 
We fit such a tail with the following form (see Fig. 1) 

P N (l) ~ A{l)U{N/N c )N^e~^ N , (35) 

where 

0(1) = \y a (i)\/N c (i), (36) 

and in general U{x) = 1. The details of the correction procedure are explained in Appendix B. 
Note that the form (35) is difficult to apply in the WR, because the CPDF naturally tends to the 
Gaussian limit as £ 2 reaches values much smaller than unity. Our method is thus practically valid 
only in the nonlinear regime, so when we correct for finite volume effects, we sample scales that 
satisfy I > Iq (see column g of Table 1). At the early stages of the simulations, where the correlation 
length is still a very small fraction of the box size, the correction for finite volume effects does not 
yield any significant change in the measured values of Sq (as expected) and is thus not necessary. 
Table 1 lists the values of a for which finite volume effects have been corrected for (column h). 

Figure 4 gives the measured values of Sq, but now corrected for finite volume effects if needed, 
in the "reliable" scaling range l m < I < ly[ (the values of the scales l m and ly[ are listed in 
Table 1). The agreement with self- similarity is considerably better than in Figure 3. The two 
n = — 1 simulations, that were giving quite different values of Sq for the uncorrected measurements 
at late stages are now in striking agreement. 

4.2.3 Measurements versus theoretical predictions 

From Figure 4, one can see a good, although not perfect, agreement of the measurements with the 
predictions of perturbation theory (long dashes) in the WR, except when n = +1, a case we further 
discuss below. Generally, the open squares on the extreme left, which correspond to the largest 
sampled scale in the first snapshot of each simulation, seem to be spurious. This is certainly due to 
finite volume effects (again) and to the fact that such weakly nonlinear scales have not yet relaxed; 
i.e. they are still influenced by transients, as discussed in § 3.2. 

In the n = +1 case, the equations (29), (30) and (31) are certainly not valid, because nonlinear 
effects may never be negligible, even when £ 2 "C 1 (this is true at least for 53, see, e.g., JBC). 
Moreover, the effective realization of the initial power-spectrum is truncated to white noise at the 
Nyquist frequency (JBC) and no power exists on scales larger than the box size, so the effective 
values of Sq in the WR should be scale-dependent functions, not necessarily close to the theoretical 
predictions (29), (30) and (31). When £ 2 1S smaller, but close to unity, the measurements are in 
agreement with the perturbation theory predictions for n = 0, as argued by JBC. For smaller values 
of £ 2 it is difficult to make any statement, since the values of Sq we measure are quite noisy, small 
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Figure 4: Same as Figure 3, but finite volume effects have been corrected for when necessary, by 
extending to infinity the large- N exponential tail exhibited by the CPDF (cf. Fig. 1); only the 
reliable scales have been displayed, i.e. l m < I < lyi, where l m and ly[ are listed in Table 1. Now, 
the agreement with self- similarity is much better: the curves corresponding to various expansion 
factors all superimpose, for a given value of Q and n. The long dashes give the predictions (29), 
(30) and (31) from perturbation theory, valid in the limit £ 2 ^ 1- I n the case n = + 1, we 
display the predictions for n = (see discussion in the text). The short dashes are the following 
phenomenological power- law fit Sq = Sq^/IOO) - 045 ^- 2 ), valid for £ 2 > 1, for all n. The values 
of Sq are given in Table 3. 

16 



or negative. Note, interestingly, that these results agree roughly with equations (29), (30) and (31) 
taken with n = +1, since in that case the predicted values of Sq are small: S3 ~ 0.86, S4 ~ 0.56 
and S 5 ~ 0.18. 

In the SR, the ratios Sq(£ 2 ) reacn a plateau, although not exactly flat. Indeed, the functions Sq 
measured in each simulation in the nonlinear regime £ 2 > 1 are a U compatible with the following 
power-law behavior 

Sq - [D(t 2 )] Q ~ 2 S Q , Q = 3,4,5, (37) 

with 

(38) 

and 

8 ~ 0.045. (39) 

The values of Sq are listed in Table 3. Note that they are quite close to those inferred from 
equations (29), (30) and (31) taken with some effective values re e fj of the power-spectral index (see 
columns e and f of Table 3). This seems to indicate that the hierarchies of correlations in the SR 
are very similar to those given by perturbation theory, at least for small Q (we neglect here the 
possible but small deviation from the scaling relation in equation [37]). The values of S3 are also in 
good agreement with the following phenomenological behavior S3 ~ 9/(3 ± n), as noticed by Fry, 
Melott & Shandarin (1993), who measured the three-body correlation function in Fourier space in 
simulations with n ranging from —3 to +1. 

The unexpected deviation from the scaling relation (5) reflected by equation (37) is very weak 
but impossible to deny, at least for n = and n = +1. Indeed, let us denote by S'g 116 the true 
values of Sq and by Sq es the measurements. The uncertainty in Sq es can be described as follows 

S t Q ™ = S Q ™f Q hl , (40) 

where fg 1 indicates that Sq es is to be multiplied by Jq or l//g. The values of Jq are listed 
in Table 3. They were estimated by considering a range of possible and reasonable values of the 
parameters involved in equation (35) and are thus valid only in the nonlinear regime £ 2 *5 1- Note 
that the fact that we have two n = — 1 realizations does not really help to make Jq closer to unity in 
that case: to have a significant improvement, we would need a few more simulations. However, we 
noticed that the measured exponential tails had quite similar features in each n = — 1 simulation, 
so we used the same parameters for both samples in equation (35) (but the very details of the 
correction are slightly different, see Appendix B). So one should not over-interpret the very good 
agreement between the two n = — 1 samples, once finite volume effects have been corrected for: 
part of this is due to our choice of the parameterization of the correction. However, if there are 
differences, they should be within /g, which corresponds to a small logarithmic shift, at most of 
order ±0.08, ±0.13, ±0.2 respectively for S3, S4 and S5. 

The dynamic range we have at our disposal in the nonlinear regime £ 2 > 1 1S a lways at least 
two orders of magnitude for the function £ 2 - Except for n = — 2, this is still the case when one 
considers the regime £ 2 <; 10. However, for the more restrictive regime £ 2 <; 100, where £ 2 agrees 
with the power-law behavior given by equations (21) and (22), we have access to such a dynamic 
range only for n = ±1. 
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Table 3: Measured values of Sq 
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a spectral index. 

b order Q of the ratio Sq. 

c the measured value of Sq when £ 2 = 100 ( see ec L- [37]). 

d the expected uncertainty factor in the measurements. If 5g ue is the true value of Sq, and 5g es 

the value we measure, then S£ ue = S^fg 1 . 
e,f the "effective" spectral index and the expected values of Sq when one tries to use formulae (29), 

(30) and (31) to fit the measured Sq in the SR. 
g The value of S3 inferred from the measurements of Fry, Melott & Shandarin (1993). 



The variations in Sq corresponding to a change of a factor 100 in £ 2 are ) according to equation 
(37) 

5 3 (100? 2 )/5 3 (? 2 ) ~ 1.2, (41) 

S 4 (100£ 2 )/S 4 (£ 2 ) - 1-5, (42) 
5 5 (100? 2 )/5 5 (? 2 ) ~ 1.85. (43) 

These ratios are close to unity but nevertheless larger than the estimated error factors Jq given 
in Table 3, except for n = — 2. There is still an ambiguity for n = — 1, but the cases n = and 
n = +1 are undoubtedly in disagreement with the scaling relation. For n = —2 and marginally for 
n = — 1, the measured Sq in the SR are compatible with the scaling relation. 

To summarize the results of this section, we find a weak but significant deviation from the stable 
clustering hypothesis for the ratios Sq(1),3 < Q < 5 in the nonlinear regime, but our measurements 
are not yet conclusive for the cases n = — 1 and n — 2. 
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5 Scaling behavior of the high-order correlations: the void prob- 
ability 



5.1 Self- similarity and spurious effects 

5.1.1 Effects of self-similarity on a 

We now want to see what the expected self-similar behavior implies for the void probability distri- 
bution function. To do so, let us rewrite equation (7) as 

a(m,l,t) = Et-l)^ 1 ^^ 11 H 3 ? 2 (W)f _1 • (44) 
N>1 

For a self-similar system, using the fact that s(t) oc a a (eqs. [14], [15]), this is equivalent to 

a(m, I, a) = a(mK~ 3 ,1k, a/c 1 /"), (45) 

where k is an arbitrary number and a is the expansion factor. To check for self- similarity, we must 
thus measure the function a(m',i',a') at number densities m! = m/c -3 ^ m. Fortunately, it is 
possible to estimate the VPDF of a realization with an average number density m' ^ m analytically 
once the function P^{m,l) is known. Indeed, we have (see, e.g., Hamilton 1985; Hamilton, Saslaw, 
& Thuan 1985; CBSII) 

oo , i\ N 

P (m,l)= E P N(n,l). (46) 

N=0 V mJ 

Note that this series converges at least for m'/m < 2. It was shown by BH and CBSII that this 
formula gives quite accurate results in practice for m'/m Sa 2; it is difficult to use for m'/m > 2. 
In the case m'/m Sa 1, the VPDF given by equation (46) will be denoted as a "diluted" VPDF. 
It is equivalent to the VPDF measured in a subsample of average number density m' randomly 
extracted from a parent sample of average number density m. 

5.1.2 Grid effects 

Figure 5 gives, in logarithmic coordinates, the quantity a as a function of l/s{t) measured in our 
simulations (except for one of the two n = — 1 samples) for various expansion factors and various 
number densities. The symbols (already used in previous figures) represent the direct measurements 
in the JV-body simulations, while the dashed, long-dashed (and dotted-dashed curves for n = —2) 
represent various diluted VPDFs using equation (46) (see the caption of Fig. 5). According to 
the self-similar solution, they should superimpose to the closest symbols, which is indeed the case, 
except perhaps for n = —2 and for the first snapshot of the n = — 1 simulation. In these cases, the 
symbols tend to lie above the dashed curves. This is mainly due to grid effects. 

Indeed, as discussed extensively in CBSII, the information carried by the regular pattern of 
particles used to set up initial conditions is likely to be conserved in underdense regions as the 
system evolves, especially if there is a cutoff at small scales in the power-spectrum (so grid effects 
should increase with — n). Consequently, the measured VPDF is smaller than it would be (so a 
is larger) if the sample was locally Poissonian. The dashed curves are not very sensitive to such 
effects because they correspond to a diluted VPDF. Indeed the prescription given by equation 
(46) is equivalent to a random dilution if m' < m: if the dilution factor k 3 is large enough, the 



19 



information linked to the regular pattern will be destroyed. The factors k 3 used in figure 5 are 
large, of order (a2/ a i) 3a ~ 16 or (a2/ a i) 6a ~ 256, where a,2 is the expansion factor corresponding 
to the parent VPDF and a\ is the expansion factor corresponding to the diluted VPDF. According 
to CBSII, the VPDF directly measured in the JV-body samples, which corresponds to the symbols 
in Figure 5 should be significantly influenced by grid effects when 

Po < 1/e- (47) 

In Figure 5, the points that satisfy Po < 1/e are circled. They indeed seem to be significantly 
further from the dashed curves for n = — 2 and for the first snapshot of the n = — 1 simulation 
(open squares). In other cases, grid effects appear to be unimportant. This is not surprising, 
because the n = +1 and n = simulations have more "power" at small scales than the n = — 2 and 
the n = — 1 simulations: larger small scale power implies earlier shell crossing and phase mixing 
which destroy the information related to the grid (see also the discussion in § 3.2). 

Note that some symbols at the extreme left of each panel also tend to lie above the dashed 
curves. We do not fully understand this small effect, but suspect it is linked to the short range 
softening of the forces. Indeed, we argued in § 4.1 that such an effect could influence £ 2 U P to 
scales of order log 10 l ~ —2.5, making the system less clustered than it should be at small scales, 
i.e. closer to a Poisson distribution for which a is unity. 



5.1.3 Finite sample effects 



The vertical segments in each panel of Figure 5 are error-bars. According to CBSII, the error in a 
can be estimated from 
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The quantity da/dN c stands for the partial derivative of <r(JV c ,m) with respect to JV C (at fixed m) 
and £2 stands for the average of £2 within the sampled volume V\> ~ r3 
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We approximate this quantity as the integral of ^2 ( | r i — r 2|) m the sphere of volume Vb OK : 

| 2 ~ ? 2 ((3/47r) 1 / 3 X box ) . (52) 



We refer to CBSII for the derivation and a discussion of the meaning and validity of equations (48), 
(49) and (50). 

There is a transition scale l cut above which AP0/P0 (and Act/it) becomes suddenly much larger 
than unity. This scale is simply defined by (CBSII) 
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Figure 5: The quantity log 10 <r(m,£) as a function of log 10 £/s, measured in each simulation for 
various expansion factors (for symbols, see Table 1) in the available scale range (eq. [11]). In 
each panel, the short-dashed curve represents an analytic dilution by a factor k 3 of the VPDF 
corresponding to the triangles using formula (46) so that, according to the self-similar solution (see 
eq. [45]), it should superimpose on the open squares. Similarly, the two long-dashed curves represent 
analytical dilutions of the VPDF corresponding to the stars. They should superimpose on the 
triangles and the open squares. The two dotted-dashed curves in the bottom right panel represent 
analytical dilutions of the VPDF corresponding to the diamonds. They should superimpose on the 
stars and the triangles. The circled points satisfy Pq < 1/e. The vertical segments are error-bars 
estimated by using equation (48). For the sake of legibility, we put the error-bars for the case 
n = —2 only on the direct measurements, i.e. the symbols. The points such that Aa/a > 1 have 
been removed. 
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Above this scale, there is typically only one "statistically independent" empty cell: if I <; l cut the 
VPDF is dominated by the largest void in the sample and its measurement is meaningless. We 
thus removed from Figure 5 the points satisfying I <; l cut . 

To estimate the value of £2 we use linear theory (eq. [20]). We can notice that the error-bars 
are insignificant for n = and n = +1, and quite reasonable for n = — 1. However, they are large 
in the n = —2 simulation, which is not surprising. Indeed, in that case (and in the case n = — 1), 
the errors are, except for the largest scales, dominated by the terms in £2 m equations (49) and 
(50). Now, the large scale "power" increases with — n, so do the values of £2 (at fixed correlation 
length). For example, for l ~ 0.13 (see Table 1) one finds 6 ^ 2 X 10" 3 , 1.5 X 10" 2 , 6 X 10" 2 and 
0.18 respectively for n = +1, 0, —1 and —2. 

Note, however, that we overestimate the finite sample error in our case from equations (48), (49) 
and (50). These equations assume that our JV-body simulations are subsamples of much bigger sets 
with fluctuations on scales larger than £t>ox- These fluctuations not only induce some uncertainty 
in the measured VPDF, but also in the average number density m through equation (49). In our 
JV-body samples, we know the average number density exactly, so the real error in our measurement 
is probably smaller than the one we compute by using the above prescription. For example, the 
difference between the values of &{m,l,a) measured in our two n = — 1 simulations are extremely 
small, much smaller than suggested by the error-bars of the bottom left panel of Figure 5 (however, 
for our example to be fully convincing, more n = — 1 realizations would be needed). Also, the 
agreement with self- similarity in the n = —2 case is excellent, well within the error-bars, except 
for the first snapshot (squares) and the largest scales, but we know that this is mainly due to grid 
effects. 

It is anyway not appropriate to set Am/m = instead of the value given by equation (49) and 
still use equation (50), because the VPDF is statistically correlated with m. By doing that, we 
would overestimate the errors even more. 

5.2 Void probability and scaling relation 

The above study of a(n,l,t) enabled us to examine self- similarity in detail to determine a dynamic 
range over which the measurements can be trusted. Now, we check whether the scaling relation 
is verified. If this were the case, a(n,l,t) would scale as a function ct(JV c ), as discussed in the 
introduction. In this system of coordinates, however, it is not easy to distinguish the WR from the 
SR, so these two regimes should be studied separately. Contrary to the regime £ 2 > 1) the range 
of values of £ 2 we probe in the regime £ 2 < 1 1S small, about one order of magnitude. Moreover, 
the measured VPDF is likely to vanish in this regime, corresponding to large cells that are likely 
to be always occupied. Simulations different from those described here would be needed to probe 
the WR. Therefore, we chose not to consider the WR in what follows. 

5.2.1 Void probability in the nonlinear regime 

Figure 6 displays the quantity a as a function of JV C measured in our simulations at various expansion 
factors a and various number densities m. The curves plotted correspond to the same values of 
(m, a) in Figure 5. We have selected the following scale range 

l>l m , Kl cut , Klo, (54) 

where l m is the "minimum reliable scale" defined in § 3 (see Table 1), and l cut is defined by equation 
(53). For the direct measurements, represented by the symbols, we moreover imposed the constraint 
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Pq > 1/e to avoid possible grid effects. This last condition is rather conservative. The results of 
§ 5.1.2 indeed indicate that grid effects on the VPDF are negligible for the cases n = and n = +1 
and are rather small in the case n = — 1, except for the smallest value of a. 

The vertical error-bars are the same as those computed in § 5.1.3. Rigorously, we should 
also have horizontal error-bars that would account for the uncertainty on the measurements of 
N c = N £ 2 - But the analysis of £ 2 m § 4.1 suggests that such uncertainties are rather small for 
I > l m , except for n = —2. In this latter case, £ 2 1S seriously affected by finite volume effects, and 
we propose a way to extract the physical information from the measurements that we use here: 
the number JV C used to construct the bottom right panel of Figure 6 is not the measured one, but 
rather the quantity 

JV c fit = JV 10 F[loglo(£/s)] , (55) 

where F is the function defined by equations (27) and (28). 

Note that for the case n = — 1, for which we have two realizations, the direct measurement of 
the quantity <r as a function of JV C gives results that are in excellent agreement for the two samples, 
except perhaps for the stars in the bottom left panel of Figure 6 that correspond to the largest 
expansion factor analyzed, a = 16. In that case, as discussed in § 4.1, £ 2 1S affected by finite volume 
effects and one can detect a small difference between the simulations. This would correspond to a 
maximum horizontal shift on the bottom left panel of Figure 6 smaller than the size of the symbols 
used to make the plots so that the measurement in the other simulation is well inside the region 
defined by the vertical error-bars. 

Clearly, for n = +1 and n = 0, a(N c ,m, a) does not scale as a function of a single variable 
ct(JV c ). The case n = — 1 is more ambiguous, but the deviation from scale-invariance is of the same 
order as the error-bars which we know overestimate the true errors. Although the measurements 
are certainly of much better quality than is suggested by the error-bars in Figure 6, the case n = —2 
is overly influenced by finite volume effects to support a definitive conclusion: the measurements 
are compatible with scale invariance as well as with a slight deviation from it. 

These results confirm those of § 4, which means that the behavior of the high order correlation 
hierarchy qualitatively follows the one of the low-order correlation hierarchy. 



5.2.2 Fine tuning 

How does the scaling behavior of the VPDF agree with the scaling behavior of the low order 
correlations? In other words, is it possible to generalize the power-law description (37) to any 
order (J? If this were the case, we would have 

&M=1 -m + m HD & )N < ) (56) 

where _ 

W)= EC- 1 )"- 1 ^'"- 1 . (57) 

N>1 

In other words, if equation (37) was valid for any Q, the quantity a = D(a — 1) + 1 would depend 
only on the number 

N c ' = D(£ 2 )N C . (58) 

This behavior is not what we find. This failure is not surprising, because equation (56) is quite 
likely to be unrealistic. Indeed, the condition Pq < 1 implies a > 0. Moreover, using a(x,y) = 
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Figure 6: Same as Figure 5, but the quantity a is represented as a function of JV C in the scaling 
range (54). Only the nonlinear scales I > Iq have been selected. For the direct measurements 
(symbols), the points satisfying Pq < 1/e have been removed (because of possible grid effects). 
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Y^N>i{~^-) N {^N{y)/N\)x N ~ 1 , we have p(N{.) = a(x = N' c ,y = l\oo), where li 00 is the scale at 
which £ 2 = 100 or equivalently D(lioo) = 1. So a must be also positive, which implies, if a < 1, 
that 

D < — ^r, (59) 

1 — (7 

a condition which is not likely to be always fulfilled, particularly when I > £100 (so D > 1) and JV C 
is large. Indeed, in the regime JV C >• 1, we find that a is always much smaller than unity in the 
reliable dynamic range (so we should have D Sa 1 in this regime). 

Thus, the tempting generalization of equation (37) to high-order correlations is not warranted. 

6 Discussion, conclusions 
6.1 Main results 

We have measured the count probability distribution function (CPDF) in a set of five high resolution 
JV-body simulations of self-gravitating expanding universes with scale-free initial power-spectra 
(|6fc| 2 ) oc k n , n = —2, —1 (two realizations), and +1. We focussed on the measurement of low- 
order averaged correlation functions £q, Q < 5 as well as high-order correlations through the void 
probability distribution function (VPDF). Our primary goal was to determine whether or not the 

— —Q_2 

expected scaling relation Sq = £q /£ 2 = constant with respect to scale was achieved in the regime 
£ 2 ^ 1) as expected if the stable clustering hypothesis applies. In our analyses, we carefully studied 
spurious effects such as grid contamination, loss of dynamics due to the short range softening of 
the forces, and sampled volume finiteness. Our main conclusions are: 

1. Finite volume effects are important for low-order correlations when Q > 3, increasingly with 
—n. Fortunately one can correct for such errors, as suggested by CBSI. We performed such 
corrections, by extrapolating to infinity the exponential tail exhibited by the CPDF at large 
N. Finite volume effects were so large in the n = —2 case in our JV-body data that they 
strongly affected the Q = 2 measurement (£2)) which was not the case for n > —1. The 
VPDF was seen to be weakly sensitive to such defects, as expected (see CBSII). 

Grid effects, due to the fact that the simulations were started from a slightly perturbed 
regular pattern of particles, are small for the low-order correlations: they affect only the 
smallest scales at early stages of the simulations, typically I ^ l c where l c is the typical 
distance between two particles in overdense regions defined by N c (l c ) = 1. Grid effects on 
the VPDF (CBSII) are significant only when n < — 1 and are of course stronger at early 
stages. They are quite strong in the case n = —2. 

In our simulations, the contamination due to the short range softening of the forces typically 
affected the measurement of the CPDF up to scales as large as 4e, where e is the softening 
length. 

2. The quantities Sq, 3 < Q < 5 considered as functions of £ 2) exhibit two plateaus separated 
by a smooth transition around £ 2 ~ 1- The measured values of Sq in the weakly nonlinear 
regime £ 2 ^ 1 (WR) are in reasonable agreement with the perturbation theory predictions. 
The values of Sq measured in the strongly nonlinear regime SR, when £ 2 > 1) are larger than 
those measured in the WR. The difference actually increases with — n. The value of Sq at 
the scale where £ 2 = 100 can be well approximated by the values predicted by perturbation 
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theory, but with an effective index re e fj = —9, —3, —1, —0.5 respectively for n = — 2, —1, 0, 
+ 1. The strongly nonlinear statistics thus behaves in a manner similar to the WR. Note also 
that the measurement of S3(£ 2 = 100) 1S m rough qualitative agreement with the findings of 
Fry, Melott & Shandarin (1993), namely that S 3 ~ 9/(3 + n). 

3. The plateau in the nonlinear regime £ 2 > 1 1S however not exactly flat, but rather of the form 

Sq oc f 2 iQ ~ 2 \ 8 ~ 0.045 for 3 < Q < 5, (60) 

with 8 ~ 0.045, independently of n. However, because of the remaining uncertainties in our 
measurements, the case n = —2 is also compatible with the scaling relation in the regime 
£ 2 *S 10. This is also true for n = — 1, but only marginally. 

4. The scaling behavior of the VPDF is in qualitative agreement with the above results: a 
significant deviation from the scaling relation is found in the nonlinear regime for the cases 
n = and n = + 1. With our error-bars, that overestimate the true errors, we see that 
our simulated n = — 1 distribution agrees only marginally with the scaling relation. The 
measurements in the n = —2 simulation are too greatly affected by finite volume effects to 
enable firm conclusion. We found that equation (60) could not be generalized to arbitrary Q. 

6.2 Discussion 

The general trend seen in our measurements is in rough agreement with the similar work of Lucchin 
et al. (1994), but with differences in detail. Indeed, our estimates are more accurate, because we 
better control various spurious effects discussed in point 1 above, particularly finite volume effects. 
Moreover, because of our high spatial resolution, we have access to a larger effective range of values 
of £ 2 , approximately one order of magnitude larger than these authors, who used a larger number 
of particles N paT = 2 192 152 than us, but evolved with a low resolution JV-body code to perform 
their simulations. However, we hardly probe the weakly nonlinear regime £ 2 Sa 1. Even if the 
comparisons of our measurements with perturbation theory are satisfying, our JV-body simulations 
are not really appropriate in that regime. More suitable investigations of the weakly nonlinear 
regime have been undertaken by Gaztahaga & Baugh (1995) for (f2, A) = (1, 0) and (0.2, 0.8) CDM 
simulations (A is the cosmological constant) and they find excellent agreement between theory and 
measurements of Sq up to Q = 10. 

The very weak deviation from the scaling relation we find in the strongly nonlinear regime 
is compatible with the earlier work of CBSI on matter distributions coming from CDM, white 
noise and HDM simulations. But it seems to disagree with some recent results in the literature, 
particularly those of Lahav et al. (1993, LIIS). If we were to believe the results of LIIS, we should 
have, in equation (60), 8 ~ 0.2 for low-density CDM and 8 ~ 0.1 for n = 0, which corresponds to 
much stronger deviations from scale-invariance than we measure. We think that the measurements 
of LIIS are contaminated by finite volume effects, particularly in the low-density CDM case. Our 
analyses, as stated by CBSI, indeed show that finite volume effects are more important when — n 
increases, or when there is a cut-off at small scales in the initial density power-spectrum, which is 
the case for a CDM spectrum when compared, for example, to a white noise spectrum. Similarly, 
Suto & Matsubara (1994) measured the three-body and four-body correlation functions £3 and £4 
(in the same simulations as LIIS) and computed various parameters associated with the hierarchical 
model (e.g., Groth & Peebles 1977; Fry 1984b) at these orders: the extrapolation of their results 
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to our averaged correlation functions suggest 8 ~ 0.1 for low-density CDM and 8 ~ 0.05 for n = 0. 
Their measurements, although of much better quality than LIIS, are certainly still contaminated 
by finite volume effects in the CDM case, but are in agreement with our analyses for white noise 
initial conditions. More recently, Bonometto et al. (1995) have measured S3 and S4 in mixed 
dark matter and in CDM JV-body samples and found strong deviation from scale-invariance. They 
attributed it to finite volume effects. Indeed, although their simulations involve a very large number 
of particles, more than 10 7 , their measurements are strongly contaminated by finite volume effects, 
mainly because the simulation box size they use (£t>ox = 100 Mpc) is not very large compared to 
the correlation length (lo -^box/10). 

To conclude this paper, as clearly illustrated by Figures 4 and 6, we establish the existence 
of a very weak but significant deviation from the scaling relation (5) in the part of the strongly 
nonlinear regime (SR) we probe, at least for n > —1. It is however important to recognize that 
our measurements do not completely reject the stable clustering hypothesis, which argues that 
the scaling relation should hold in the regime £ 2 ^ 1 if a self-similar solution is reached. Indeed, 
earlier measurements by EFWD show that the power-law behavior (21) of £ 2 implied by the stable 
clustering hypothesis is achieved only when £ 2 *S 100, as illustrated by Figure 2. In this regime, we 
hardly probe an order of magnitude of values of £ 2 f° r n = —1 and n = —2, a bit more for n = 
and two orders of magnitude for n = +1. With this more restrictive approach, the deviation from 
the scaling relation we measure is really significant only in the case n = +1, and it is anyway quite 
small. One would need much larger high resolution simulations to significantly probe the strongly 
nonlinear regime £ 2 , with large number of particles (say more than 128 3 instead of our 64 3 ) and 
higher spatial resolution (smaller short range softening length e). 
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APPENDIX 

A Timesteps 

Following Efstathiou et al. (1985, EDFW), we use the following time variable: 

p = a a , a=— ^— , (61) 
F n+3 y ' 

where a is the expansion factor, taken to be unity at the beginning of the simulation. In the units 
of the treecode of Hernquist et al. (1991), the gravitational constant is G = 1 and the mass of a 
particle is m = N~^ T . To probe the same dynamic range as EFWD, all the simulations were evolved 
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from p = 1 to p max = 16, except for n = —2 where p max = 64. The corresponding values of the 
expansion factor are given in Table 1. 

According to EDFW, the timestep Ap should satisfy the following requirement 

1/2 

(62) 

at any time for the time-centered leapfrog scheme to be stable against round-off errors, where rj 
is the softening parameter used by EDFW in their P 3 M code. This rj corresponds to an effective 
softening length (EFWD) 

e ~ 77/3 (63) 

for a potential of the form 

</> oc (r 2 + e 2 )~ 1/2 . (64) 

Although, we use a slightly different softening of the short range forces (i.e. a cubic spline interpo- 
lation; see Hernquist & Katz 1989; Goodman & Hernquist 1991), we assume that eq. (63) is still 
valid in our case. Equation (62) can now be rewritten 

Ap <S 1.73aa a (e/A par ) 3 / 2 , (65) 

where A par is the mean interparticle distance. 

For e = A par /20, which is the value we choose for the softening parameter, we find 

Ap & 0.0194aa a , (66) 

which leads to Ap & 0.0097, 0.013, 0.019, 0.039 for n = 1,0,-1,-2 respectively. Such val- 
ues also apply to the scale invariant simulations of EFWD and are indeed in very good agree- 
ment with the values of Ap displayed on Table 1 of EFWD. Our choice of timestep was Ap = 
0.0084, 0.011,0.017,0.034 respectively for n = 1,0, -1, -2, which is slightly smaller than the one 
imposed by the condition (66). 



Ap <S ^atLa a+1 ' 2 



Gm 



B Correcting for finite volume effects 

This paragraph explains in detail how we use the phenomenological form (35) to correct for finite 
volume effects. 

Since £ 2 1S n °t affected by finite volume effects, except for the case n = —2 which requires 
a special treatment as explained below, we use the measured values of JV C in equation (35). We 
replace the measured CPDF by the analytical form (35) for (arbitrary) N larger than a lower 
bound iV cr ; t (^). This number iV cr ; t must be large enough so that the function indeed exhibits 
the behavior implied by equation (35). Typically, we should have iV cr ; t a few JV C when U(x) = 1 
(but practically, we take JV crit £ N c ) and JV crit > N™£ with JV c ^ t n > 1 (e.g., CBSI). The last 
condition is imposed because of discreteness effects. In some cases for n = +1 and n = —2, we used 
a more elaborate fit than equation (35) with U = 1, to have a description as accurate as possible 
of the shape of the CPDF. In that case, iV cr ; t can be smaller than N c . Following BSD and BH, we 
take 

U(x) = (1 + bx)~ c , (67) 
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Table 4: The values of the parameters used to correct for finite volume effects 
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40 



lo S10 I 
^crit l^c 

^crit /*c 

^crit /"c 
_V 

^crit l^c 

^crit /^c 
J7 

^crit /^c 
_>7 



1.5 
-2.6 

1.0 
-2.45 
0.05 
-1.7 

2.0 
-2.1 

1.5 
-2.1 

2.0 
-2.1 



1.5 
-2.6 

1.0 
-2.45 
0.05 
-1.7 

2.0 
-2.2 

1.5 
-2.1 

1.5 
-2.1 



1.5 
-2.6 

1.0 
-2.45 
0.05 
-1.7 

2.0 
-2.1 

1.5 
-2.1 

1.3 
-2.1 



1.5 
-2.6 

1.0 
-2.45 
0.05 
-1.7 

2.0 
-2.2 

1.5 
-2.1 

1.1 
-2.1 



2.0 
-2.8 
1.0 
-2.45 
0.05 
-1.7 

2.0 
-2.35 
1.5 
-2.1 
0.98 
-2.1 



2.5 
-3.0 

1.0 
-2.45 
0.05 
-1.7 

2.2 
-2.7 

1.5 
-2.2 

0.7 
-2.2 



1.0 
-2.57 
0.05 
-1.7 



1.0 
-2.68 



1.5 
-2.7 

1.9 
-2.7 



-"crit 
V 



/JVc 



2.5 
-1.575 
2.0 



2.5 
-1.575 
2.0 



2.5 
-1.575 
2.0 



2.5 
-1.575 
2.0 



2.3 
-1.575 
1.8 



3.0 
-1.575 
1.5 



-"crit 
V 



/JVc 



3.0 
-1.575 
2.2 



2.0 
-1.575 
2.0 



2.0 
-1.575 
2.0 



2.0 
-1.575 
2.0 



1.0 
-1.575 
2.0 



1.0 
-1.575 
2.0 



1.0 
.575 



2.0 
-1.575 
1.5 



1.7 
-1.35 



1.7 
-1.35 



1.7 
-1.35 



1.7 
-1.35 



1.7 
-1.35 



1.7 
-1.35 



1.7 
-1.35 



.7 
.35 



1.7 
-1.35 



a spectral index. 
b expansion factor. 

c the correction parameters are the same for both n = — 1 simulations at a = 6.4. 
d the value of N CT - lt /N c is slightly different for each of the n = — 1 runs at a = 16. 
e in the case n = 0, it appears appropriate to have scale dependent \y s \ and constant rj. 



where b and c are a priori adjustable parameters that may vary with scale, but that can be fixed 

at constant values at the level of precision adopted here. To follow CBSI, since £ 2 1S n °t affected 

by finite volume effects, except for n = —2, we impose a choice of the parameters such that the 

correction does not change the value of £ 2 - 

The case n = —2 is slightly more complicated. For a = 3.2, where finite volume effects are still 

small for the function £ 2 ) we use the above procedure. For larger values of a, we use the fitting 

— _ At 

function F(x) determined in § 4.1.3 to describe £ 2 an d then write N c = N £ 2 m equation (35) with 

£fit = 1Q F[log 10 (//.)]. (68) 

We require that the measured function £ 2 , once corrected for finite volume effects, be as close as 
possible to £ 2 . 

Table 4 lists the values of iV cr ; t /iV c , iV™ t n , r](l), \y s (l)\, b, c we chose when we corrected for finite 
volume effects. Note that when 6 = and c = 1, this simply means we took U(x) = 1. 
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